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We examine the electronic properties of 2D electron gas in black phosphorus multilayers in the 
presence of a perpendicular magnetic field, highlighting the role of in-plane anisotropy on various 
experimental quantities such as ac magneto-conductivity, screening, and magneto-plasmons. We 
find that resonant structures in the ac conductivity exhibits a red-shift with increasing doping due 
to inter-band coupling, 7. This arises from an extra correction term in the Landau energy spectrum 
proportional to n 2 y 2 (n is Landau index), up to second order in 7. We found also that Coulomb 
interaction leads to highly anisotropic magneto-excitons. 

PACS numbers: 72.80.Vp,85.85.+j,73.63.-b 


I. INTRODUCTION 


The successful exfoliation of black phosphorus (BP) 
multilayers has triggered tremendous interests in this 
material, which is also one of the thermodynamically 
more stable phases of phosphorus, at ambient temper¬ 
ature and pressure. In comparison to other 2D materi¬ 
als, such as graphene, hexagonal boron nitride (hBN), 
and transition metal dichalcogenides, multilayers BP has 
a direct bandgap which spans from 0.3 eV ~ 1.5 eV 5 ®, 
hence making it an excellent candidate for infrared 
optoelectronics 7 "®. Since each BP layer forms a puckered 
surface due to sp 3 hybridizat ion, i t also reveals highly 
anisotropic electrical mobility 1 3 3®, linear dichroism in 
optical absorption spectrePI 7 23G2] anis otrop ic excitonic 
structure 1 ®^] an d anisotro pic p lasmond 1 ®®®. However, 
BP is not stable in ambient 17 ^® which might render its 
electrical properties less than pristine. 

Recently, encapsulation of BP with hexagonal boron 
nitride (hBN)P®, all within a controlled inert atmo¬ 
sphere, has allowed for higher carrier mobility in these 
BP devices 1 ®? ^Sl. Similar hBN encapsulation has also 
been applied to other 2D materials such as graphene®^ 
and transition metal dichalcogenideiP® to achieve record 
mobilities. Indeed, the high quality BP has made 
possible the first observation of prominent quantum 
magneto-oscillations in these devices 19 23 and quantum 
Hall effect 7 . Hence, theoretical studies of BP in the 
presence of a magnetic field has also begun receiving 

attentiorPEB 7 . 

In this paper, we examine the electronic properties 
of BP multilayer thin film in perpendicular magnetic 
field, such as its Landau level spectrum, ac conductiv¬ 
ity, screening and its collective electronic excitations. In 
particular, we emphasize the manifestation of anisotropy 
in these experimentally observable quantities. We begin 
with a discussion of the model Hamiltonian used to de¬ 
scribe multilayer BP in Section II. This is followed by 


the study of its electronic subband structure in Section 
III, and its Landau level spectrum in Section IV. Various 
experimentally relevant quantities such as ac conductiv¬ 
ity, collective excitations and screening in the presence of 
magnetic field will be discussed in Section V-VII respec¬ 
tively. 


II. MODEL HAMILTONIAN 


In multilayer BP, broken translational symmetry in the 
out-of-plane z direction renders the direct energy gap at 
the T point instead of the Z point in bulk. The low 
energy Hamiltonian description of BP near T point can 
be expressed as TL = TL Z + TL xy , with its out-of-plane and 
in-plane dynamics taken separately. Here, TL XV is given 
bjBra 
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xy 
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where r] c ^ v and v c , v are the respective band parameters, 
while 7 describe the effective coupling between the con¬ 
duction and valence bands. E c and E v denotes the ener¬ 
gies of the bulk conduction and valence band edges. We 
discuss our choice of these band parameters below. 

Cyclotron resonance experiments on bulk BF 29 found 
an out-of-plane electron and hole effective masses con¬ 
siderably smaller than that of layered tansition metal 
dichalcogenides materialiP®. In this work, we adopt an 
average of experimental 2 -® and theoretically^^® predicted 
out-of-plane masses i.e. m cz ss 0.2mo and m vz « 0.4mo, 
mo being the electron mass. The out-of-plane Hamilto¬ 
nian is given by, 


_fi 2 / m c 3 d 2 
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eV(z) 


( 2 ) 


where V(z) describes the out-of-plane electrostat ic po- 
tential, typically induced by a bottom metal gat#®!*®®. 
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FIG. 1: Electronic subband structure of 10 nm BP film for 
electron (a) and hole (b), showing the lowest four subbands, 
5Ec' V . The red dashed lines corresponds to the Fermi energy 
Ef. These calculations assume T = 300 K. 



For a finite BP thickness with given V(z), Eq. ([2| can be 
diagonalized numerically, leading to electronic subband 
structure. We denote these subband eigen-energies as 
5E 3 V , and their eigen-functions as <jP cv {z), where j is the 
subband index. 

Close to the T point, the in-plane band parameters, 
i.e. 7) c ,v, v c v and 7 , are related to the in-plane effective 
masses vi J®, 


h 2 

= ±2 ( 7 2 /E g + r 1(c , v) ) 


c,v)y 


h 2 


(3) 


The band parameters rj c>v , u c/v and 7 are chosen such 
that they yield the known effective masses in the bulk 
BP limit i.e. m cx = m vx — 0.08 too, m cy = 0.7 too and 
m vy = 1.0 to.( J I 0 I 29 ( and m cx = m vx ss 0.15 Too for mono- 
layer BFGSI. Eg is the electronic bandgap of the BP mul¬ 
tilayer. The energy gap for bulk BP is 0.3 eV 1 ^ While 
monolayer BP has not been ascertain experimentally, ab 
initio calculation based on the GW method suggests an 
energy gap of ~ 1.5 — 2e’W^. In sum, the band pa¬ 
rameters used in this work are; i] c = r] v ss 19.1 eV A 2 , 
v c « 5.45eVA 2 , u v « 3.81 eVA 2 and 7 = 2.84eVA. 


III. SUBBAND STRUCTURE 

Typical experimental device structure 19 ^ consists of 
multilayer BP on an insulating dielectric film on a back 
substrate, which serves also as a back gate, as sketched 
in the inset of Fig.jlja). The electron density in BP along 
2 , n(z), can be obtained froirPQ, 


kgT ^ 

n{z) = smcd x 


In 


exp 


E f - E c - SEi 
kT 
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where g s is the spin degeneracy, kg is the Boltzmann 
constant, Ef is the Fermi level, T is the temperature, 


FIG. 2: Landau level spectrum for the conduction (a) and 
valence (b) subbands, calculated for a 10 nm thick BP film. 
The underlying subband structure corresponds to the case 
with a carrier density of 1 x 10 12 cm -2 . In (c), we show the 
probability distribution ip(£ = s/ax , \fo.y) for the two low¬ 
est Landau levels in conduction band under different Landau 
gauges at B = 10T, where a = oj c /2hri c . In (d), the proba¬ 
bility distribution of the i — 0, 5, 10, 15 eigenstates over the 
bare oscillator eigenstates | n) of lowest conduction subband 
is plotted. 


and m c d refers to the density-of-states mass given by 
sjm cx m cy . Solving Eq. (|2j) and the Poisson equation self- 
consistently, one can then arrives at the numerical solu¬ 
tion for the BP electrostatics, an approach well-known in 
the context of semiconductor inversion layei^l. 

The electron and hole subbands, 8E 3 CV , and the Fermi 
energy Ef , are plotted as function of carrier density n 
in Fig.[l](a)-(b) respectively. These calculation assumes 
T = 300 K. The results indicate that across a wide range 
of carrier densities, only the first subband is occupied, 
and onset of second subband occupation takes place only 
when n > 7 x 10 12 cm -2 and n > 5 x 10 12 cm -2 for elec¬ 
tron and hole respectively. This is consistent with recent 
experimental observation&A21®3 of a 2D electron gas in 
the quantum limit. Certainly, the transition to multi¬ 
subband occupation depends also on the BP thickness. 


IV. LANDAU LEVEL SPECTRUM 


When an uniform magnetic field Bz is applied perpen¬ 
dicular to the plane, we have p = Hk —> II = p + eA by 
Peierls substitution. With the choice of Landau gauge 
A = Bxy , we can write 

n x - Pxi Lly - Py ± CXB) 


(5) 
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which obey the commutation relation [ 11 ^, n y ] = — iHeB. 
It is useful to define, 


a = 


1 

y/flUJ c 


[y/Vc^-x 


1y/ ^clly] 


( 6 ) 


where ui c = 2eB^/rj c v c . It is easy to verify that [a, a^] = 
1 , and 

at ° = E7^ cI1 ^ + Vc ^ ^ 

Introducing the parameters s = r/ v /g c , t = u v /u c and 
r = s/t, we can express the dispersing term in valence 
band as: 


where 7 = f - r ) and e iv = 

(E' cv + 5E 0 c v )/huj c . We note that, for the representa¬ 
tive numerical results presented here, the band edges 
E c v are adjusted so as to reproduce the estimated elec¬ 
tronic bandgap of ss 0.5 eV^ for a 10 nm BP him, i.e. 
{E' c + 8El)~{E' v +5El) = 0.5 eV. 

Now the matrix Hamiltonian is dimensionless, and the 
eigenvalue problem can be solved numerically. The pro¬ 
cedure is as follows. First, we assume the following ansatz 
for the eigenvector of Eq. ([9]) , 



En =0 U n (j)n{cH{x X 0 )) 
E n=0 X Vn0n{a{x - X 0 )) 


.( 10 ) 


ri v U 2 x + v v Il 2 y = s(r] c Il 2 x + p c n 2 y ) + (1 - r)v v n 2 (8) 

which allows us to obtain the effective Hamiltonian un¬ 
der magnetic held B for each pair of electron-hole j sub¬ 
bands, 

v j = h,, ( e c + at ° + § . ^ + a ) ^ 

\ 7( a ^ + a ) ej + s{a i a + |) + r(a) — a) 2 J 


where n max sets the truncation of the expansion and 
4> n {a{x — xo)) is the n’th eigenstate of harmonic oscil¬ 
lator centered at x 0 = with a = . Note the 

eB 2nr] c 

wave function will include e lPvV ^ h factor, associated to 
the good quantum number p y . We can determine the 
coefficients U n and V n from the energy eigen-problem 
(9) Wip = Eif). Explicitly, we have the following relations 


(Ai + n — E)U n + jiy/riVn—i + yjn + lV^+i) — 0 
^{y/nUn -1 + y/n + lU n+ i) + (A 2 + sn - E)V n + f(y/n(n- 1 ) 14-2 + V( n + 2 )( n + 1)^+ 2 ) = 0 


( 11 ) 


where Ai = A 2 = + | — r, s = s — 2r and 

E = E/huj c , E being the energy of Landau level cor¬ 
responding to the eigen-function (10). Numerically, we 


need to introduce a truncation condition, which we set 

to U n> n mlzn: = Vn>ro maa . = 0 * 

It is not difficult to check that for most physical rel¬ 
evant cases, A 1)2 n,j, which justifies a perturba¬ 
tive consideration for Eq. (11). For the electron (hole) 


spectrum, we can write, E e / h = E e J h + SE e ^ h with 
E 0 e = Ar + n, E% = A 2 + sn, and the second order 
perturbations terms are: 


5E e / h = ±7 2 (c 0 + c\n + c 2 n 2 ) 


( 12 ) 


where cq = 


and c 2 


! 

(Ai — A 2 ) ’ 
~ 2(1 - 5) 
(Ar — A 2 ) 2 


1 — s 


Cl = 


(Ai — A 2 ) (Ai — A 2 )^ 


From Eq. (12), we may under¬ 


stand the deviation of resonant frequency of the ac con¬ 
ductivity from conventional 2D electron gas case, to be 
discussed in Section V. 

From numerical recipe described below Eq. ( |11|), we ob¬ 
tain the Landau level spectrum, as shown in Fig.r2[a)-(b) 
for conduction and valence bands respectively. The dis¬ 
persion is typical of 2D electron gas system, exhibiting 


linear dependence with B. However, there are differences 
due to the finite inter-band coupling 7 and anisotropy of 
BP, which will be discussed later. In the calculations, 
we set the doping to be 1 x 10 12 cm -2 . We note that 
crossing of each Landau level acquires additional carrier 
density of An = eBg/h ~ \B x 10 11 cm -2 , where g = 2 
is the spin degeneracy. Hence, at the assumed doping 
of 1 x 10 12 cm -2 and B = 10 T, the filling factor is 2. 
For carrier densities larger than 1 x 10 13 cm' 2 , the filling 
factor can be as large as 20. It might then be possible 
to observe multi-subband phenomena, especially for the 
hole case. For the experimental doping range, only the 
Landau levels of the lowest two subbands are physically 
relevant. 

The anisotropy of the problem is encoded in the wave- 
functions. In Fig.[2](c) we plot the wavefunction probabil¬ 
ity along the two in-plane spatial coordinates for the first 
two Landau levels, expressed in their dimensionless coor¬ 
dinate (i.e., £ = y/ax, £ = y/ay). Two respective gauges 
are used, i.e., for gauges A = Bxy and A = —Byx , we 
get ip(x) and V'(y) respectively. Due to the anisotropy 
inherent in the model Hamiltonian, we can clearly dis¬ 
cern the difference between probability distribution over 
these two gauges. We will show that this anisotropy in 
wave function can result in prominent anisotropy in vari- 
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ous experimental quantities such as ac conductivities and 
magneto-plasmons. In Fig.0d), the probability distribu¬ 
tion over the eigenstates | n) of the bare harmonic oscil¬ 
lator which forms the upper diagonal terms in Eq. (|9 ) is 
plotted. With the increase of Landau level, the probabil¬ 
ity distribution over |n) becomes more broadened. This 
can be understood from perturbation point of view, i.e., 
the matrix elements quantifying the perturbation upon 
the bare eigenstates increases with factors y/n’s. 

Before concluding this section, it is interesting to com¬ 
pare Landau levels in BP with the other electron gas sys¬ 
tem e.g. conventional 2D electron gas (i.e. Schrodinger 
fermions) and graphene (i.e. Dirac fermions), summa¬ 
rized as follows: 

{ hoj c (n + |) : Schrodinger fermions 
sgn(n)vfy/2ehB\n\ : Dirac fermions (13) 
huj c {n + \ + 7 2 (c 0 + Cin + c 2 n 2 )) : BP 

The BP effective Hamiltonian, described by Eq. 0 , have 
features not embodied in Schrodinger fermion description 
in the finer energy scale proportional to q 2 , and have ex¬ 
perimental consequences discussed in Section V. We note 
also that compared to other gapped Dirac systems, such 
as gaped graphene or transition metal dichalcogenides, 
BP’s gap is placed at the time-reversal invariant T point 
instead of the inequivalent K and K' points of the BZ. 
Hence, its Landau level spectrum resembles more to that 
of Schrodinger fermions rather than that of massive Dirac 
fermions . 7 


150 - 





(a) B= 10 T 



100 - 





—r 1 - 

- n = 2 >?,'/ V, o 

- 0->1 


—V 


- 1-»2 

50 - 


n= 0 

- 2->3 





0 - 

■ ■ m m 4 ■-& 





tl Co (meV) 


FIG. 3: ac conductivities as function of frequency to. The 
anisotropy between a xx and a yy , at B = 10 T is displayed 
in (a), for different n —> n + 1 transitions. Zooming in, we 
can observe slight shift of the resonant frequency for a xx and 
a xy , as shown in (b), calculated for B = 20 T. The damping 
constant is set to be V = 10~ 4 eV and temperature T = 10 K. 
The inset of (b) shows how the resonant frequency shift de¬ 
pends on filling factor v and B. Conductivity are in units of 
cto = e 2 /h. 


V. AC MAGNETO-CONDUCTIVITY 

In this section, we study the ac magneto-conductivity 
a a p of the 10 nrn BP multilayer him. Having the en¬ 
ergy spectrum and wave-functions of Landau levels at 
hand, as presented in the previous section, we can nu¬ 
merically calculate the cr a p’s directly. According to the 
Kubo formulfP^, we have: 


Ca/sM = 


X 


-ie 2 h f{E sn j ) - f{E a i n iji) 

iii E sn j E s f n fjt 
snjs'n'j' J J 

(^*snj |^q| ^s'n'j') s'rij ' \vf}\ ‘F snj) 

Hto E sn j -P E s i n 'j' T fr 


where the velocities are defined by Vi = ® r -, Explicitly, 
we have: 


v x = 


( y]hw c rj c (a + a 1 ') 7 

7 Vvy/f^c/Vcia + a 1 ') 




( y/hio c v c (a) — a)/i 

l 0 


(15) 


,y/hu c /v c (at - a)/i J 


Our results for ac conductivity as a function of fre¬ 
quency are presented in Fig. [3] The ac longitudinal 


magneto-conductivities for filling factor v = 1 to 3 (elec¬ 
tron doped) is shown in Fig.[3](a). The inset depicts the 
transitions between the nearest Landau levels for cases 
with v = 1 to 3, i.e. the resonance for each case corre¬ 
sponds to a particular transition process |n) to | n + 1 ), 
and occurs at the terahertz frequencies. Prominently, 
a xx is about 5 ~ 10 times larger than <j yy . We note the 
spatial anisotropy in the wavefunctions of the Landau 
levels, albeit small, as shown in Fig. 0 - The anisotropy 
in the magneto-optical conductivity tensor arises mainly 
from the anisotropy in the velocity operators, i.e. v x and 
v y , which accounts for the anisotropic optical transitions 
dipole. 

We should point out here that the resonant structure 
observed here in BP is more like conventional 2D elec¬ 
tron gas rather than in graphene^, which has multiple 
resonant structures. Contrary to the conventional case, 
we find here that the resonant frequency is slightly red- 
shifted with increasing doping (or v). This is shown in 
Fig.§b), where the longitudinal conductivity a xx and 
Hall conductivity a xy are displayed for io = 1—4, cal¬ 
culated for B = 20 T. This red-shift also increases with 
magnetic field as depicted in the inset. We note that 
the frequency shift increase in almost uniform steps each 
time the filling factor decreases by 1, Interestingly, this 
red-shifting behavior can be understood from the per- 
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FIG. 4: Excitation spectrum of BP, as obtained by density 
plots of ImII(q, oj). (a) and (b) correspond to non-interacting 
polarization of BP, Eq. ( |17[ ), for q || y and q || x respec¬ 
tively. Plots (c) and (d) include electron-electron interactions 
in the RPA. The strength of the electron-electron interaction 
is chosen to be r 3 ~ 3, and Nf = 3. 


turbation expression in Eq. (12 1 rather straightforwardly, 
from which we arrived at, 


E e {n + 1) - E e (n) = 1 + y 2 [ci + c 2 (2n + 1)] (16) 


This expression accounts for linear dependence of the res¬ 
onance frequency on n quantitatively, and is a direct re¬ 
sult of the interband coupling 7 . In other words, this 
red-shift in ac conductivity can be used to determine 7 
in BP. As depicted in the inset of Fig.[3](b), it can be 
seen that the linear relation holds when v is small. We 
anticipate experimental progress in magnetic oscillations 
of BP to shed light on this issu^lMl ft j s a l so worthy 
to mention that the longitudinal and Hall conductivities 
calculated here can be directly measured through absorp¬ 
tion and Faraday rotation experiments 34 via terahertz 
spectroscopy. 


VI. COLLECTIVE EXCITATIONS 

In this section we study the excitation spectrum of 
BP in the presence of a quantizing magnetic field ap¬ 
plied perpendicular to the sample, including the effect 
of electron-electron interaction. As we have seen in Sec. 
|III[ a 2DEG is created in a multi-layer BP, with occupa¬ 
tion of only the first subband, unless the doping exceeds 
~ 5x 10 12 cm -2 . Therefore in this section we concentrate 


on the particle-hole excitation spectrum of the 2DEG 
formed by the carriers of a lOnm thick BP multi-layer. 
In the absence of Coulomb interaction, the particle-hole 
excitation spectrum, that enclose the region of the uj — q 
plane in which it is possible to excite electron-hole pairs, 
can be calculated from Imn°(q, uj) 7 ^ 0 , where n°(q, uj) 
is the non-interacting polarization function. The bare 
polarizability of BP in the quantum Hall regime can be 
expressed in terms of the standard result for a 2DEGP 


n°(q,w) = Y1Y! 


n,m (q) 

ui — mu) c + ir 


+ (w + -» —w ) (17) 


where = E„=max(o ,N F -m) and w+ indicates 

the replacement ui + *r —» — u> — zT, and Np is the in¬ 
dex of the last occupied LL. As we have seen before, 
the anisotropy of the BP band structure is encoded in 
the wave-function, whose overlaps leads to different form 
factors in the x- and y-direction 


E n,m(q) — 


n! 


27 xl 2 B 


(n + to)! 


a- 


2;2 

B 


qH 


n 2 I 2 
q ij. 


( 18 ) 

where l b = \Jhc/eB is the magnetic length, a = 
Lo c /2hrj( c v ) for the case when q x is a good quantum num¬ 
ber, and a = uj c /2hv( c ^ when the gauge chosen leads 
to a good q y quantum number. We note that the finite 
interband coupling, 7 , can renormaliz^H the band pa¬ 
rameters (? 7 ( Ci „)and V( c , v )) but we checked that the effect 
is small in this case. A density plot of Imn° is shown in 
Fig. [4|a) and (b) for electron doping with Np = 3. In 
the presence of a quantizing magnetic field, Imn°(q,w) 
is a sum of Lorentzian peaks centered at u> = mu> Cl where 
m is the difference between the LL indices of the electron 
n! and the hole n, m = n! — n > l ! 36 * 37 The excitation 
spectrum is chopped into horizontal lines, separated by a 
constant energy u c . The peculiarities of the BP spectrum 
of Fig. |4](a)-(b), like the presence of a superstructure of 
Np + 1 brighter regions and the nodes of the first hor¬ 
izontal line at uj = w c , are due to the form of the LL 
wavefunctions and have been studied in detail in Refs. 
1371551 As in the B = 0 casej^ anisotropy of BP band 
structure leads to a wider spectrum in the q y direction as 
compared to q x direction. Whereas we only show results 
for electron doping in Fig. |4j a similar spectrum is found 
for the hole doped case, with the difference that the sep¬ 
aration between consecutive horizontal lines is narrower 
for the latter, due to larger effective mass of the valence 
band as compared to the conduction band. 

We next consider the effect of Coulomb interaction 
within the Random Phase Approximation (RPA), which 
dresses the electron-hole polarization function as 


n(q,w) = 


n°(g^) 

e(q,w) 


(19) 


where the dielectric function is obtained as 


e(q,w) = 1 - u(q)n°(q,w), 


( 20 ) 
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in terms of the non-interacting polarization function II 0 
and the two-dimensional Coulomb potential in momen¬ 
tum space 

, . 27re 2 

*>(q) = —j—r ( 21 

«|q| 

where k is the background dielectric constant. The 
strength of the Coulomb interaction is usually ex¬ 
pressed in terms of the dimensionless parameter r s = 
2mbe 2 /nkp, where the band mass in the present case is 
mb = yjm x m yi and \zf = \j2Np + 1/Ib = y/Ttn e i in 
terms of the carrier density n e j. Long-range electron- 
electron interaction leads to the appearance of collective 
plasmon modes in the spectrum. Their dispersion rela¬ 
tion can be calculated from the zeros of the dielectric 
function. In Fig. |4jc) and (d) we show how the non¬ 
interacting electron-hole spectrum [shown in Fig. EM 
an (b)] is modified due to interactions. The electron-hole 
horizontal lines of Fig. [4j[a) and (b) acquire a dispersion, 
leading to the so-called magneto-excitons or magneto- 
plasmons. As in the B = 0 caseP^ Coulomb interac¬ 
tion leads to highly anisotropic magneto-plasmons, with 
higher dispersion for q || x than for q || y directions. 
This is due to the fact that the effective mass is smaller 
in the x direction than in the y direction. We also no¬ 
tice that the excitation spectrum of BP greatly differs 
from that of doped graphene. In fact, the characteris¬ 
tic linear dispersion relation in graphene leads to a rel¬ 
ativistic quantization of the spectrum into a set of non- 
equidistant LLsP^ As a consequence, long range Coulomb 
interaction in graphene leads to a set of highly dispersing 
linear magneto-plasmonspSl that differs from the equidis¬ 
tant magnetoexcitons (separated by a well defined cy¬ 
clotron frequency energy, ca c ) in BP. 

VII. STATIC SCREENING 

In this section we focus on the properties of II 0 (q) = 
II 0 (q, uj = 0) and e(q) = e(q,w = 0) in the static 
limit, for which the polarization function is entirely real. 
The polarizability of BP in a strong magnetic field is 
shown in Fig. 51a). One first observe that, as in a stan¬ 
dard 2DEGp3"the static polarizability tends to zero as 
II 0 (q -A 0) oc q 2 for B ^ 0. The reason for this is that 
the main contribution to II 0 (q) comes from q = 0 excita¬ 
tions in the vicinity of the Fermi energy Ef. This differs 
from the B = 0 case, for which the Fermi level cuts the 
band and there are q = 0 excitations. In the integer 
quantum Hall regime, however, Ef lies in the cyclotron 
gap between the highest occupied LL Np and the low¬ 
est unoccupied one Np + 1. Since this energy gap must 
be overcome by excitations with q = 0 , then its spectral 
weight tends to zero. Notice that n°(q = 0) coincides 
with the density of states at the Fermi energy because 
the latter vanishes for B > 0 when Ef lies in the gap. 
One also notices that the wave-vector at which the po¬ 
larizability starts to vanish (2k p) is larger for q |j y that 


-n°[q,o] 




FIG. 5: Static screening. Static polarization function 
—n°(q, 0 ) (top) and dielectric function e(q) (bottom) for the 
conduction band with q || y (solid lines) and for q || x (dashed 
lines). We have used Np = 3. 


for q || x, due to the anisotropy of the BP Fermi surface. 
One further notices the oscillatory behavior of the static 
polarizability, below 2k p, due to the wave-function over¬ 
lap bet ween t he electron and the hole, leading to Np + 1 
maxima ! 36 1 37 1 

It is interesting to compare the screening properties of 
BP and graphene, and we start by briefly discussing the 
B = 0 case. For q — > 0 we have e(q) «1 + 5 tf/ 9 , where 
qxF = 2ire 2 p(Ef)/n is the 2D Thomas-Fermi wave- 
vector, in terms of the density of states p(Ef). Since the 
density of states (per unit area) is approximately a con¬ 
stant equal to g s mb/2n for BP, where g s = 2 accounts for 
the spin degeneracy, whereas for graphene p(Ef) is en¬ 
ergy dependent and given approximately by gEf / (2i TVp), 
where vp is the Fermi velocity and g = g s g v = 4 accounts 
for spin and valley degeneracy, one obtains a density inde¬ 
pendent qpF for BP, whereas it scales as kp for graphene. 
Therefore, in the two cases and for B = 0, the dielectric 
function diverges as e ~ qrF/q -A oo when q -A 0. 15 
However, it is important to notice that the density de¬ 
pendence in the numerator of e(q —» 0) in graphene im¬ 
plies the absence of screening in undoped graphene (i.e. 
for kp = 0) at long distances. For short wavelengths 
q 2 kp, e(q) —> 1 in a BP, whereas for graphene, 
e(q) -A 1 + irr s /2. The extra contribution 7 rr s /2 to the 
dielectric function of graphene at q 3> 2k f is related 
to virtual inter-band particle-hole excitations.^ In sum¬ 
mary, at short wavelengths and for B = 0, BP (like a 
standard 2DEG) does not screen at all (e —► 1), whereas 
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graphene screens as a dielectric due to its filled valence 
band. 

The above picture changes in the presence of a quan¬ 
tizing magnetic field. In Fig. [5] we have plotted the 
static polarization and dielectric functions for BP, re¬ 
spectively. At long wavelengths, e(q) — 1 oc q. In the 
limit of Np 1, the known result for a 2DEG and for 
graphene,^ s(q) — 1 oc r s N^/ 2 qlB as q —¥ 0, applies also 
for BP. The difference between graphene and BP is en¬ 
coded in the density dependence of r s in the two cases: 
since r s ~ N F ' in BP, e grows linearly with Np in 
this case. However, r s = e 2 /kvf is density-independent 
in graphene, leading to a dielectric function proportional 
to Ny 2 . On the other hand, in both graphene and BP, 
e(q,uj = 0) —► 1 as q —»• 0, which implies that there is 
no screening at long distances. Finally, the behavior of 
the dielectric function in a magnetic field at g > 2 kp in 
BP is similar to the B = 0 limit, corresponding to the 
standard metallic like screening governed by intra-band 
processes.^ 

VIII. CONCLUSION 

In conclusion, we have examined the electronic prop¬ 
erties of 2D electron gas in black phosphorus multilayers 


due to the presence of a perpendicular magnetic field. We 
highlight in this work the in-plane anisotropy reflected in 
various experimental quantities such as its ac magneto¬ 
conductivity, screening, and collective excitations. We 
found that resonant structures in the ac conductivity ex¬ 
hibits a red-shift with increasing doping due to interband 
coupling, suggesting possible electric modulation of light 
absorption and Faraday rotation. Coulomb interaction 
also leads to highly anisotropic magneto-plasmons. 
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